

In [1]:

    
%matplotlib inline
from pylab import *
from numpy import *



In [2]:

    
#physical constants
hbar = 1.0
m = 1.0



In [3]:

    
#discretization
nx = 100

minx = -1.
maxx = 1.

x = linspace(minx, maxx, nx)



In [4]:

    
#assemble the matrix Hamiltonian
from d2matrix import d2matrix
from scipy.sparse import dia_matrix

def make_H(V):
    nx = len(V)
    kinetic = -hbar**2/(2*m)*d2matrix(nx)
    potential = diag(V)
    H = kinetic + potential
    return H


V = zeros(nx)

H = make_H(V)









    



/usr/local/lib/python3.4/dist-packages/numpy/core/fromnumeric.py:2507: VisibleDeprecationWarning: `rank` is deprecated; use the `ndim` attribute or function instead. To find the rank of a matrix see `numpy.linalg.matrix_rank`.
  VisibleDeprecationWarning)



In [5]:

    
# calculate the eigenvectors
eigvals, eigvecs = linalg.eig(H)
eigvecs = array(eigvecs)



In [6]:

    
for i in arange(nx):
    plot(x, eigvecs[:,i]+eigvals[i]*150)
ylim([0, 10])









    Out[6]:





(0, 10)



In [7]:

    
def solve(H):
    eigvals, eigvecs = linalg.eig(H)
    
    I = eigvals.argsort()
    eigvals = eigvals[I]
    eigvecs = array(eigvecs)[:,I]
    return eigvals, eigvecs



In [8]:

    
V = x**2
H = make_H(V)
eigvals, eigvecs = solve(H)



In [9]:

    
figure(figsize(5,6))
for i in arange(10):
    plot(x, eigvecs[:,i]*.05+eigvals[i])
plot(x, V, lw=3)
ylim([0, 0.3])









    Out[9]:





(0, 0.3)

Two atoms:



In [10]:

    
#discretization
nx = 100

minx = -1.
maxx = 1.

x = linspace(minx, maxx, nx)

V = zeros(nx)
a1 = -0.5
a2 = 0.5
asize = 0.2
apot = -0.07
V[(x>a1-asize) & (x<a1+asize) | (x>a2-asize) & (x<a2+asize) ] = apot



In [11]:

    
eigvals, eigvecs = solve(make_H(V))



In [12]:

    
figure(figsize(6,6))
for i in arange(6):
    plot(x, eigvecs[:,i]*.05+eigvals[i])
plot(x, V, lw=2)
set_printoptions(precision=16)
print(eigvals[:6])

#ylim([-1, -0.6])









    



[-0.0624387366052539 -0.062438578859082  -0.0406528132169814
 -0.0406493136442565 -0.0091801384671852 -0.0089414977975879]

Evolution of a localized electron? Evolution of each eigenvector follows from Schrodinger equation $$ i\hbar\frac{\partial\psi_j}{\partial t} = \hat H\psi_j = E_j\psi_j $$ $$ \psi_j(t) = e^{{E_j t}/{i \hbar}}\psi_j(0) $$ And for general initial state $$ \psi(0) = \sum_j c_j\psi_j(0), {\rm where }\quad c_j = \langle\psi_j(0)|\psi(0)\rangle $$ is given by $$ \psi(t) = \sum_j c_j\psi_j(t) = \sum_j c_j e^{{E_j t}/{i \hbar}}\psi_j(0) $$



In [15]:

    
coeffs = zeros(nx)
coeffs[0] = 1/sqrt(2.)
coeffs[1] = 1/sqrt(2.)
E_mean = dot(coeffs*eigvals, coeffs)



In [16]:

    
def psi_t(eigvecs, eigvals, coeffs, t):
    phase = exp(eigvals*t/1j*hbar)
    #print phase
    
    psi = zeros(len(coeffs), dtype=complex)
    for i in arange(len(coeffs)):
        if coeffs[i] != 0.:
            psi += eigvecs[:,i]*phase[i]
    return psi



In [17]:

    
from matplotlib import animation
from JSAnimation.IPython_display import display_animation



In [18]:

    
fig = figure()
ax = axes(xlim=(minx,maxx))
line1, = ax.plot([],[],lw=1)
ax.plot(x, V)


def animate(data):
    line1.set_data(x, abs(psi_t(eigvecs, eigvals, coeffs, 2000000.*data))*.1+E_mean)

    return line1,

anim = animation.FuncAnimation(fig, animate, frames=arange(20), interval=100)
display_animation(anim, default_mode='loop')

http://www.hep.manchester.ac.uk/u/forshaw/BoseFermi/Double%20Well.html



In [ ]: